function [] = plot_disp_map(TQ_tab,plotmogi,plotmt)
%josh crozier 2020

%DEM
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %import dem
% DEM = struct();
% [DEM.Z,DEM.x,DEM.y] = import_hawaii_dem();
% domain = [257,265,2143.5,2151.5]*1000;
% DEM.Z = DEM.Z(:,DEM.x>=domain(1)&DEM.x<=domain(2));
% DEM.x = DEM.x(DEM.x>=domain(1)&DEM.x<=domain(2));
% DEM.Z = DEM.Z(DEM.y>=domain(3)&DEM.y<=domain(4),:);
% DEM.y = DEM.y(DEM.y>=domain(3)&DEM.y<=domain(4));
% 
% %downsample
% DEM.Z = imgaussfilt(DEM.Z,2);
% DEM.downsamp = floor(length(DEM.x)/100);
% DEM.x = DEM.x(1:DEM.downsamp:end);
% DEM.y = DEM.y(1:DEM.downsamp:end);
% DEM.Z = DEM.Z(1:DEM.downsamp:end,1:DEM.downsamp:end);


%station data
VariableNames = TQ_tab.Properties.VariableNames;
stationcell = {};
for ii = 1:length(VariableNames)
    if length(VariableNames{ii}) > 7
        if strcmp(VariableNames{ii}(end-3:end), '_HHZ')
            if strcmp(VariableNames{ii}(1:4), 'amp_')
                stationcell{end+1} = VariableNames{ii}(8:10);
            end
        end
    end
end
stations = struct('station',stationcell);

%station locations
for ii = 1:length(stations)
    if strcmp(stations(ii).station,'NPB')
        stations(ii).latlon = [19.41200 -155.28100 1115.0];
    end
    if strcmp(stations(ii).station,'SRM')
        stations(ii).latlon = [19.39531 -155.27400 1128.0];
    end
    if strcmp(stations(ii).station,'NPT')
        stations(ii).latlon = [19.41200 -155.28100 1115.0];
    end
    if strcmp(stations(ii).station,'OBL')
        stations(ii).latlon = [19.41750 -155.28400 1107.0];
    end
    if strcmp(stations(ii).station,'WRM')
        stations(ii).latlon = [19.40650 -155.30000 1163.0];
    end
    if strcmp(stations(ii).station,'SDH')
        stations(ii).latlon = [19.38950 -155.29400 1133.0];
    end
    if strcmp(stations(ii).station,'UWE')
        stations(ii).latlon = [19.42111 -155.29300 1240.0];
    end
    if strcmp(stations(ii).station,'UWB')
        stations(ii).latlon = [19.42469 -155.27800 1091.0];
    end
    if strcmp(stations(ii).station,'SBL')
        stations(ii).latlon = [19.42700 -155.26800 1084.0];
    end
    if strcmp(stations(ii).station,'HAT')
        stations(ii).latlon = [19.42300 -155.26100 1082.0];
    end
    if strcmp(stations(ii).station,'BYL')
        stations(ii).latlon = [19.41200 -155.26000 1059.0];
    end
    if strcmp(stations(ii).station,'KKO')
        stations(ii).latlon = [19.39781 -155.26600 1146.0];
    end
    if strcmp(stations(ii).station,'RIM')
        stations(ii).latlon = [19.39531 -155.27400 1128.0];
        stations(ii).station = 'RIMD';
    end
    if strcmp(stations(ii).station,'PUH')
        stations(ii).latlon = [19.385510 -155.251310 1079];
        stations(ii).station = 'PUHI';
    end
end
%convert lat-lon locations to UTM E-N coordinates
zone = utmzone(stations(ii).latlon(1),stations(ii).latlon(2));
utmstruct = defaultm('utm');
utmstruct.zone = zone;
utmstruct.geoid = wgs84Ellipsoid;
utmstruct = defaultm(utmstruct);
for ii = 1:length(stations)
   [stations(ii).locE, stations(ii).locN, stations(ii).locZ] = ...
       mfwdtran(utmstruct,stations(ii).latlon(1),stations(ii).latlon(2),stations(ii).latlon(3)); 
end

%lava lake location (from Chouet work)
% epilatlon = [19.405402 -155.281041 100.0];
% [epiloc_east, epiloc_north, ~] = mfwdtran(utmstruct,epilatlon(1),epilatlon(2),epilatlon(3));

if plotmogi
    %mogi location (from Kyle Anderson Geodetic inversions for 2018 eruptions)
    mogi_east = 63.7 + 260740;
%     eastoffset = kyle_east-epiloc_east;
    mogi_north = -41.95 + 2147480;
%     northoffset = kyle_north-epiloc_north;
end

if plotmt
    %mogi location (from Kyle Anderson Geodetic inversions for 2018 eruptions)
    mt_east = 63.7 + 260740;
%     eastoffset = kyle_east-epiloc_east;
    mt_north = -41.95 + 2147480;
%     northoffset = kyle_north-epiloc_north;
end

%displacements amp,phase
for ii = 1:length(stations)
    in = find(strcmp(strcat('amp_HV_',stations(ii).station,'_HHZ'), VariableNames));
    stations(ii).amp_HHZ = TQ_tab{1,in};
    in = find(strcmp(strcat('amp_HV_',stations(ii).station,'_HHE'), VariableNames));
    stations(ii).amp_HHE = TQ_tab{1,in};
    in = find(strcmp(strcat('amp_HV_',stations(ii).station,'_HHN'), VariableNames));
    stations(ii).amp_HHN = TQ_tab{1,in};
    in = find(strcmp(strcat('phase_HV_',stations(ii).station,'_HHZ'), VariableNames));
    stations(ii).phase_HHZ = TQ_tab{1,in};
    in = find(strcmp(strcat('phase_HV_',stations(ii).station,'_HHE'), VariableNames));
    stations(ii).phase_HHE = TQ_tab{1,in};
    in = find(strcmp(strcat('phase_HV_',stations(ii).station,'_HHN'), VariableNames));
    stations(ii).phase_HHN = TQ_tab{1,in};
    if plotmogi
        in = find(strcmp(strcat('amp_HV_',stations(ii).station,'_HHZ_mogidepth'), VariableNames));
        stations(ii).amp_HHZmogi = TQ_tab{1,in};
        in = find(strcmp(strcat('amp_HV_',stations(ii).station,'_HHE_mogidepth'), VariableNames));
        stations(ii).amp_HHEmogi = TQ_tab{1,in};
        in = find(strcmp(strcat('amp_HV_',stations(ii).station,'_HHN_mogidepth'), VariableNames));
        stations(ii).amp_HHNmogi = TQ_tab{1,in};
        in = find(strcmp(strcat('phase_HV_',stations(ii).station,'_HHZ_mogidepth'), VariableNames));
        stations(ii).phase_HHZmogi = TQ_tab{1,in};
        in = find(strcmp(strcat('phase_HV_',stations(ii).station,'_HHE_mogidepth'), VariableNames));
        stations(ii).phase_HHEmogi = TQ_tab{1,in};
        in = find(strcmp(strcat('phase_HV_',stations(ii).station,'_HHN_mogidepth'), VariableNames));
        stations(ii).phase_HHNmogi = TQ_tab{1,in};
    end
    if plotmt
        in = find(strcmp(strcat('amp_HV_',stations(ii).station,'_HHZ_mtdepth'), VariableNames));
        stations(ii).amp_HHZmt = TQ_tab{1,in};
        in = find(strcmp(strcat('amp_HV_',stations(ii).station,'_HHE_mtdepth'), VariableNames));
        stations(ii).amp_HHEmt = TQ_tab{1,in};
        in = find(strcmp(strcat('amp_HV_',stations(ii).station,'_HHN_mtdepth'), VariableNames));
        stations(ii).amp_HHNmt = TQ_tab{1,in};
        in = find(strcmp(strcat('phase_HV_',stations(ii).station,'_HHZ_mtdepth'), VariableNames));
        stations(ii).phase_HHZmt = TQ_tab{1,in};
        in = find(strcmp(strcat('phase_HV_',stations(ii).station,'_HHE_mtdepth'), VariableNames));
        stations(ii).phase_HHEmt = TQ_tab{1,in};
        in = find(strcmp(strcat('phase_HV_',stations(ii).station,'_HHN_mtdepth'), VariableNames));
        stations(ii).phase_HHNmt = TQ_tab{1,in};
    end
end

keepin = isfinite([stations.amp_HHZ] + [stations.amp_HHE] + [stations.amp_HHN]);
stations = stations(keepin);
VariableNames = VariableNames(keepin);

%reference phase
ref_in = find(strcmp({stations.station},'NPT'));
if isempty(ref_in)
    ref_in = find(strcmp(stations.station,'NPB'));
end
if isempty(ref_in)
    ref_in = 1;
end
ref_phase = stations(ref_in).phase_HHZ;
for ii = 1:length(stations)
    stations(ii).phase_HHZ = stations(ii).phase_HHZ - ref_phase;
    stations(ii).phase_HHE = stations(ii).phase_HHE - ref_phase;
    stations(ii).phase_HHN = stations(ii).phase_HHN - ref_phase;
    if plotmogi
        stations(ii).phase_HHZmogi = stations(ii).phase_HHZmogi - ref_phase;
        stations(ii).phase_HHEmogi = stations(ii).phase_HHEmogi - ref_phase;
        stations(ii).phase_HHNmogi = stations(ii).phase_HHNmogi - ref_phase;
    end
    if plotmt
        stations(ii).phase_HHZmt = stations(ii).phase_HHZmt - ref_phase;
        stations(ii).phase_HHEmt = stations(ii).phase_HHEmt - ref_phase;
        stations(ii).phase_HHNmt = stations(ii).phase_HHNmt - ref_phase;
    end
end

%displacement vectors
t = linspace(0,2*pi,80);
for ii = 1:length(stations)
    stations(ii).Z = stations(ii).amp_HHZ*cos(stations(ii).phase_HHZ);
    stations(ii).E = stations(ii).amp_HHE*cos(stations(ii).phase_HHE);
    stations(ii).N = stations(ii).amp_HHN*cos(stations(ii).phase_HHN);
    if plotmogi
        stations(ii).Zmogi = stations(ii).amp_HHZmogi*cos(stations(ii).phase_HHZmogi);
        stations(ii).Emogi = stations(ii).amp_HHEmogi*cos(stations(ii).phase_HHEmogi);
        stations(ii).Nmogi = stations(ii).amp_HHNmogi*cos(stations(ii).phase_HHNmogi);
    end
    if plotmt
        stations(ii).Zmt = stations(ii).amp_HHZmt*cos(stations(ii).phase_HHZmt);
        stations(ii).Emt = stations(ii).amp_HHEmt*cos(stations(ii).phase_HHEmt);
        stations(ii).Nmt = stations(ii).amp_HHNmt*cos(stations(ii).phase_HHNmt);
    end
    stations(ii).Zlist = stations(ii).amp_HHZ*cos(stations(ii).phase_HHZ + t);
    stations(ii).Elist = stations(ii).amp_HHE*cos(stations(ii).phase_HHE + t);
    stations(ii).Nlist = stations(ii).amp_HHN*cos(stations(ii).phase_HHN + t);
    stations(ii).Ecirc = stations(ii).amp_HHZ*cos(t);
    stations(ii).Ncirc = stations(ii).amp_HHZ*sin(t);
    if plotmogi
        stations(ii).Zlistmogi = stations(ii).amp_HHZmogi*cos(stations(ii).phase_HHZmogi + t);
        stations(ii).Elistmogi = stations(ii).amp_HHEmogi*cos(stations(ii).phase_HHEmogi + t);
        stations(ii).Nlistmogi = stations(ii).amp_HHNmogi*cos(stations(ii).phase_HHNmogi + t);
        stations(ii).Ecircmogi = stations(ii).amp_HHZmogi*cos(t);
        stations(ii).Ncircmogi = stations(ii).amp_HHZmogi*sin(t);
    end
    if plotmt
        stations(ii).Zlistmt = stations(ii).amp_HHZmt*cos(stations(ii).phase_HHZmt + t);
        stations(ii).Elistmt = stations(ii).amp_HHEmt*cos(stations(ii).phase_HHEmt + t);
        stations(ii).Nlistmt = stations(ii).amp_HHNmt*cos(stations(ii).phase_HHNmt + t);
        stations(ii).Ecircmt = stations(ii).amp_HHZmt*cos(t);
        stations(ii).Ncircmt = stations(ii).amp_HHZmt*sin(t);
    end
end

%plot
figure()
% contour((DEM.x)/1000,(DEM.y)/1000,DEM.Z,0:20:3000,'LineColor',[0.5,0.5,0.5])
% hold on
scale = 1.5/max(sqrt([stations.E].^2 + [stations.N].^2));
scatterscale = 20;
for ii = 1:length(stations)
%     plot(stations(ii).locE/1000 + stations(ii).Elist*scale,...
%         stations(ii).locN/1000 + stations(ii).Nlist*scale,'b')
    p1 = plot(stations(ii).locE/1000 + stations(ii).Ecirc*scale,...
        stations(ii).locN/1000 + stations(ii).Ncirc*scale,'b');
    hold on
%     plot(stations(ii).locE/1000 + [0,0],...
%         stations(ii).locN/1000 + [0,stations(ii).amp_HHZ]*scale,'-c+')
    if plotmogi
%         plot(stations(ii).locE/1000 + stations(ii).Elistmogi*scale,...
%             stations(ii).locN/1000 + stations(ii).Nlistmogi*scale,'r')
        p1mogi = plot(stations(ii).locE/1000 + stations(ii).Ecircmogi*scale,...
            stations(ii).locN/1000 + stations(ii).Ncircmogi*scale,'r');
%         plot(stations(ii).locE/1000 + [0,0],...
%             stations(ii).locN/1000 + [0,stations(ii).amp_HHZmogi]*scale,'-m+')
    end
    if plotmt
%         plot(stations(ii).locE/1000 + stations(ii).Elistmt*scale,...
%             stations(ii).locN/1000 + stations(ii).Nlistmt*scale,'r')
        p1mt = plot(stations(ii).locE/1000 + stations(ii).Ecircmt*scale,...
            stations(ii).locN/1000 + stations(ii).Ncircmt*scale,'g');
%         plot(stations(ii).locE/1000 + [0,0],...
%             stations(ii).locN/1000 + [0,stations(ii).amp_HHZmogi]*scale,'-m+')
    end
end
p2 = quiver([stations.locE]/1000,[stations.locN]/1000,...
    [stations.E]*scale,[stations.N]*scale,'b','AutoScale','off','LineWidth',2);
% quiver([stations.locE]/1000,[stations.locN]/1000,...
%     zeros(1,length(stations)),[stations.Z]*scale,'c','AutoScale','off','LineWidth',2)
if plotmogi
    p2mogi = quiver([stations.locE]/1000,[stations.locN]/1000,...
        [stations.Emogi]*scale,[stations.Nmogi]*scale,'r','AutoScale','off','LineWidth',1);
%     quiver([stations.locE]/1000,[stations.locN]/1000,...
%         zeros(1,length(stations)),[stations.Zmogi]*scale,'m','AutoScale','off','LineWidth',1)
    psourcemogi = plot(mogi_east/1000,mogi_north/1000,'rs','MarkerSize',9);
end
if plotmt
    p2mt = quiver([stations.locE]/1000,[stations.locN]/1000,...
        [stations.Emt]*scale,[stations.Nmt]*scale,'g','AutoScale','off','LineWidth',1);
%     quiver([stations.locE]/1000,[stations.locN]/1000,...
%         zeros(1,length(stations)),[stations.Zmt]*scale,'m','AutoScale','off','LineWidth',1)
    psourcemt = plot(mt_east/1000,mt_north/1000,'gs','MarkerSize',9);
end
if plotmogi && ~plotmt
    legend([p1,p2,p1mogi,p2mogi,psourcemogi],["VLP Vertical","VLP Horizontal",...
        "Mogi Vertical","Mogi Horizontal","Mogi Source"],'AutoUpdate','off','Location','NorthWest')
elseif ~plotmogi && plotmt
    legend([p1,p2,p1mt,p2mt,psourcemt],["VLP Vertical","VLP Horizontal",...
        "MT Vertical","MT Horizontal","MT Source"],'AutoUpdate','off','Location','NorthWest')
elseif plotmogi && plotmt
    legend([p1,p2,p1mogi,p2mogi,psourcemogi,p1mt,p2mt,psourcemt],["VLP Vertical","VLP Horizontal",...
        "Mogi Vertical","Mogi Horizontal","Mogi Source",...
        "MT Vertical","MT Horizontal","MT Source"],'AutoUpdate','off','Location','NorthWest')
end
plot([stations.locE]/1000,[stations.locN]/1000,'k.','Markersize',9)
text([stations.locE]/1000 + 0.1,[stations.locN]/1000,{stations.station})
xlabel('East (km)')
ylabel('North (km)')
grid on
box on
axis equal
xlim([257.5,263.5])
ylim([2144.5,2151])
set(gca,'ydir','normal')
end

